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Abstract 



Using computational algebraic geometry techniques and Hilbert bases of polyhedral cones we derive explicit 
formulas and generating functions for the number of magic squares and magic cubes. 

Magic cubes and squares are very popular combinatorial objects (see |l^, ^ and their references). 
A magic square is a square matrix whose entries are nonnegative integers and whose row sums, column 
sums, and main diagonal sums add up to the same integer number s. We will call s the magic sum of 
the square. In the literature there have been many variations on the definition of magic squares. For 
example, one popular variation of our definition adds the restriction of using the integers 1, . . . as 
entries (such magic squares are commonly called natural or pure and a large part of the literature consists 
of procedures for constructing such examples, see ||, ^), but in this article the entries of the squares 
will be arbitrary nonnegative integers. We will consider other kinds of restrictions instead: 

Semi-magic squares is the case when only the row and column sums are considered. This apparent 
simplification has in fact a very rich theory and several open questions remain (see Q Q and references 
within. Semi-magic squares are called magic squares in these references). Pandiagonal magic squares 
are magic squares with the additional property that any broken-line diagonal sum adds up to the same 
integer (see Figure ||) . 
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Figure 1 : Four broken diagonals of a square and a pandiagonal magic square. 



There are analogous definitions in higher dimensions. A semi-magic hypercube is a d-dimensional 
n X n X ■ ■ ■ X n array of n"^ non-negative integers, which sum up to the same number s for any line 
parallel to some axis. A magic hypercube is a semi-magic cube that has the additional property that the 
sums of all the main diagonals, the 2'^~-^ copies of the diagonal xi^i^...^!, a;2,2,....2, • ■ ■ ,Xn,n,...,n under the 
symmetries of the d-cube, are also equal to the magic sum. For example, in a 2 x 2 x 2 cube there are 4 
diagonals with sums xi^i^i + X2,2,2 — X2,i,i + 2^1,2,2 — 2;i,i,2 + a;2,2,i — a;i,2,i + ^2,1,2- We can see a magic 
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Figure 2: A magic cube. 



Two fundamental problems about magic arrays are (1) enumerating such arrays and (2) generating 
particular elements. In this paper we address these two issues from a discrete-geometric perspective. 
The work of Ehrhart and Stanley |l^, ^ |2^ when applied to the study of semi-magic squares 
showed that many enumerative and structural properties of magic arrays can actually be formulated in 
terms of polyhedral cones. The conditions of constant magic sum can be written in terms of a system 
{a;|Aa: = 0, a; > 0}, where the vector x has as many entries as there are cells in the array (labeled 
Xi-^^i2,...,ia)i E^nd a matrix A with entries 0, 1 or —1 forces the different possible sums to be equal. 

The purpose of this note is to study the convex polyhedral cones defined by magic squares, pandiag- 
onal magic squares, semi-magic hypercubes, and magic hypercubes. In particular we study the Hilbert 
bases and extreme rays of these cones. We have used computational polyhedral geometry and commu- 
tative algebra techniques to derive explicit counting formulas for the four families of magic arrays we 
defined. Similar derivations had been done earlier for semi-magic squares |^ , §4] . The interested reader 
can download the complete extre me ray information and Hilbert bases from [www . math . ucdavis . edii/ 
"deloera/RESEARCH/magic . html| 

Hilbert bases for these cones of magic arrays are special finite sets of nonnegative integer arrays that 
generate every other nonnegative integer array as a linear nonnegative integer combination of them. Most 
of our arguments will actually use minimal Hilbert bases which are smallest possible and unique ||2^ . 
Due to their size and complexity, our calculations of Hilbert bases and extreme rays were done with the 
help of a computer. We explain later on our algorithmic methods. 

Having a Hilbert basis allows the generation of any magic array in the family, and makes trivial the 
construction of unlimited numbers of such objects or simply to list all magic arrays of fixed small size. 
Another benefit is that a Hilbert basis can be used to compute generating functions for the number of 
magic arrays from the computation of Hilbert series of the associated affine semigroup ring. We carry on 
these calculations using Grobner bases methods. Finally minimal integer vectors along extreme rays of a 
cone are in fact also members of the Hilbert basis. 

It is well-known from the work of Ehrhart that for any rational pointed cone, if its lattice points 
receive a grading (e.g. by total sum of the entries, or in this case magic sum), then the function that counts 
the lattice points of fixed graded value is a quasipolynomial. A function / : N ^ C is quasipolynomial if 
there is an integer iV > and polynomials /o, . . . , /n-i such that f{m) = fi{m) if m = i (mod N). The 
integer TV, which is not unique, will be called a quasi-period of /. If it is the smallest quasi-period it will be 
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called the period of the quasipolynomial. It is a natural question to investigate when the quasipolynomial 
is actually a polynomial, i.e when the period is one. We study this question for the four families of 
magic arrays. It is known that for the cone of semi-magic squares the quasi-polynomial is actually a 
polynomial (i.e. period is one). This follows from a well-known result of Ehrhart that assures that, for 
integral polytopes, the function that counts lattice points inside their integral dilations is a polynomial. 
We prove that the same argument does not work for other magic arrays. This will involve studying the 
extreme rays of the various cones. In general, determining exactly the period is a delicate issue as seen 
from Example 4.6.27 [0. 

Consider the convex hull P of real nonnegative arrays (of given size) all whose mandated sums equal 
1. We call the polytope P the polytope of stochastic magic arrays. For example, the stochastic semi-magic 
squares are the well-known bistochastic matrices (n x n matrices whose row and column sums are one) 
and P is the famous Birkhoff-von Neumann polytope |^, It is easy to see that the polytope P can 
be written as P = {x G i?'' : a; > 0, and Bx = 1} where the matrix B has {0, 1} entries. B has as many 
rows as axial sums (row, column, diagonals, etc), and the columns of B correspond to the entries of the 
magic array. 

In Bona presented a proof that the counting function of semi-magic 3x3x3 cubes is a quasi- 
polynomial of non-trivial period. In our first theorem we extend this by actually computing an explicit 
generating function and quasipolynomial formulas for the number of semi-magic 3x3x3 cubes. 

Theorem 0.1. Denote by SHf^{s) the number of semi-magic d- dimensional hypercubes with n'^ entries. 
We have the following results 

1. From the Hilbert bases for the cones of 3 x 3 x 3 semi-magic cubes we obtain the generating function. 
EZo SHl{s)t' = *'+5i^+67t'^ + 130«^_+242«^+130«^^^^^^ 1 + 12t + 132^^ + SATt^ + 3921<4 + 

14286^5 + 43687t6 -f 116757i^ -I- . . . ). 
In other words, 



{ 9 I <;7 I =6 _L 297 „5 , 1341 „4 , 513 „3 , 3653 „2 , 627 „ , i -j- 91 « 

2240 ^ 560 "^ 320 320 ^ ^ 640 160 ^ 1120 ^ 280 ^ '■J ^P' 

9 ,8 -I- ^ -I- JL <j6 I 297 5 I 1341 =4 , 513 3 , 3653 2 , 4071 „ i nfhprim-ip 

2240 ^ 560 ^ 320 * ^ 320 ^ 640 * ^ 160 ^ 1120 * ^ 2240 ^ 128 ^ * 

2. The number of vertices of the polytope of stochastic semi-magic n x n x n cubes is bounded below by 
(n!)^"/?!" . The polytopes of stochastic 3x3x3x3 semi-magic 3x3x3 cubes and 3x3x3x3 
hypercubes are not integral. 

We also computed an explicit generating function for the number of 3 x 3 x 3 magic cubes. 

Theorem 0.2. Let AfC„(s) denote the number of n x n x n magic cubes. Then, AfC„(s) is a quasipoly- 
nomial of degree {n — 1)^ — 4 for n > 3, n ^ 4. For n — A it has degree (4 — 1)'^ — 3 = 24. For n = 3, 
using the minimal Hilbert basis for the cones o/ 3 x 3 x 3 magic cubes, we computed X]^o -^^'3(5)^^ = 
ti^+i4ts'+36t«+i4t^+i I + I9t^ I21t^ + A39t^ + mit'^'^ 2581t^^ -\- A999t^^ + . . .). Thus, in terms 
of a quasipolynomial formula we have: 
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r Jis4 , n 3 , 25 2 I 7 , -, 
MCsis) = <^ 324* ^ 54* ^ 36* 6 * + ^ ""I*' 

otherwise. 
The polytope of stochastic 3x3x3x3 magic hypercubes is not integral. 

Our next contribution is to continue the enumerative analysis done in [Q. These authors wrote down 
formulas for the number of magic squares of orders 3 and 4. We have corrected a minor mistake in the 
4x4 formula of page 8] (the 3x3 case has been known since 1915 ^^), we find further values for 
order 5 magic squares and we give evidence supporting one of their conjectures page 9]. 

Theorem 0.3. // Af„(s) denotes the number of n x n magic squares of magic sum s, then , from the 
minimal Hilbert bases for the cones o/ 4 x 4 and 5x5 magic squares, we obtain 

EZoM4{s)t' = *'+"*'+^^*°+ffi'+4fi*!+.^f'+^'*'+^*+^ (= l + 8i + 48i2 + 200t3 + 675t4 + 1904t5+4736t6 + 
10608i^ + 21925i« + ...;, 
specifically we obtain that 



480 * ^ 240 ^ 480 * ^16 ^ 30 * ^15 ^ 30 * ^ ^1*' 

J_s7,_^ 6, ^ 5, 11 4, 779 3, 593 2, 1051, 13 otherwise 

480 ^ 240 ^ 480 * ^16 ^ 480 ^ 240 * ^ 480 * ^ 16 * 

We also know the values of Mc,{s) for s < 6. The polytope of stochastic magic squares is not integral 
for 71 > 2. 

Finally, we continue the work started in ^ for the study of pandiagonal magic squares. Here 
we investigate their Hilbert bases, as an application we recomputed the formulas of Halleck (see |]l6| . 
Chapters 8,10]. The integrality of the polytope of panstochastic magic squares was fully solved in 

Theorem 0.4. Let AfP„(s) denote the number ofnxn pandiagonal magic squares of magic sum s, then 
from the Hilbert bases for the cones of 4 x 4 and 5x5 pandiagonal magic squares we obtain 



MP4(.)^/ ^(.s2+4. + 12)(. + 2)2 zf2\s, 



) otherwise. 

^'^-^5(5) = + 4)(s + 3)(s + 2)(s + l)(s2 + 5s + 8)is'' + 5s + 42). 

8064 

Here is the plan for the paper: In Section ^ we review the notion of (minimal) Hilbert bases and how 
we computed them. We show how to use a Hilbert basis to compute a generating function that counts 
the number of nonnegative integer arrays of given magic sum. In that section we recall some basic facts 
about polyhedral cones, Ehrhart polynomials, and commutative semigroup rings (see p4] , [2^). Finally, 
in Section |^, we discuss the specific details for the four theorems above, each appearing in a separate 
subsection. We close this introduction remarking that the algebraic-geometric techniques used here are 
not the only useful computational tools. In fact, there has been a surge of interest on such techniques 
with good practical results (see [sj ^). 
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1 Hilbert bases for counting and element generation 



Let A be an integer d x n matrix, we study pointed cones of the form C — {x\Ax = 0,x > 0}. A cone is 
pointed, if it does not contain any linear subspace besides the origin. It is weU-known that pointed cones 
admit also a representation as the set of all possible nonnegative real linear combinations of finitely many 
vectors, the so called extreme rays of the cone (see page 232 of (24|). As an example we consider the cone 
of 3 X 3 magic matrices. This cone is defined by the system of equations 



Xll + Xi2 + Xi3 = X21 + X22 + X23 = X31 + X32 + X33 

xii + xi2 + xi3 = Xll + X21 + X31 = xi2 + X22 + X32 a;i3 + X23 + X33 

Xll + X12 + Xi3 = Xll + X22 + X33 = X31 + X22 + Xi3, 



and the inequalities xij > 0. In our example for 3x3 magic squares the cone C has dimension 3, it is 
a cone based on a quadrilateral, thus it has 4 rays (see Figure It is easy to see that all other cones 
that we will treat for magic arrays are also solutions of a system Ax = 0, a; > 0, where A is a matrix with 
0, 1, — 1 entries. For a given cone C we are interested in Sc — C fl Z", the semigroup of the cone C. 

An element v of Sc is called irreducible if a decomposition v = vi + V2 for vi,V2 G Sc implies that 
= or V2 = 0. A Hilbert basis for C is a finite set of vectors HB{C) in Sc such that every other 
element of Sc is a positive integer combination of elements in HB{C). A minimal Hilbert basis HB[C) 
is inclusion minimal with respect to all other Hilbert bases of C. As a consequence all elements of the 
minimal Hilbert basis HB{C) are irreducible and HB{C) is unique. 
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Figure 3: The Hilbert basis for the cone of 3 x 3 magic squares. The top four squares are the rays of the 
cone. 



A natural question is then, how can we compute the minimal Hilbert basis of a cone C ? Several 
research communities have developed algorithms for computing Hilbert bases having different applications 
in mind: integer programming and optimization p^ , commutative algebra |Q, and constraint 

programming [ pO| , p^ . In our calculations of minimal Hilbert bases we used extensively the novel project- 
and-lift algorithm presented in and implemented in MLP by R. Hemmecke. On the other hand we were 
able to corroborate independently most of our results using a different algorithm, the cone decomposition 
algorithm, implemented in NDRMALIZ by Bruns and Koch 0. Similar ideas were also discussed in ||27|] . 
Now we present brief descriptions of these two methods. 
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Hemmecke's algorithm for computing the Hilbert basis of a pointed rational cone C expressed as 
{z : Az = 0, z £ M" } proceeds as follows: Let Hj : M" W be the projection onto the first j coordinates. 

Let Kj -.^ {7r]'(w) : v e kerz(A)}, K+ := Kj n {R^^ x and KJ^ := Kj n {R^^ x R-), where 

kerz{A) — {z : Az — 0,z ^ Z"} denotes the integral kernel (or null space) of A. Observe that Kj~ 
and K~ are semi-groups under vector addition. Let Hj' and H~ denote the unique inclusion minimal 
generating sets of the semi-groups {K^ ,+) and {K^^, +). Clearly, H = since K:^ = C. 

The idea of the project-and-lift algorithm is to start with , which is easy to compute, and to 
compute iJjt^i U H~^^ from Hj~ . This last step is done by a completion procedure (similar to s-pair 
reduction in Buchberger's algorithm and is based on the fact that for any vector v S kerz{A) with 
TTj^i{v) G U the vector TTj{v) can be written as a non-negative integer linear combination of 

elements in Hj' . Since many unnecessary vectors are already thrown away when Hj'^^ is extracted from 
Hj]^^ U H~^i, intermediate results are kept comparably small and larger problems can be solved. 

The cone-decomposition algorithm, used in NORMALIZ triangulates the cone C into finitely simplicial 
cones. A cone is simplicial if it is spanned by exactly n linearly independent vectors ui, . . . , w„. There 
are many possible triangulations, and any of these can be used. For each simplicial cone consider the 
parallepiped 11 = {XiVi -I- • • ■ -I- A„u„ G Z"|Ai G [0, 1), }. It is easy to see that the finite set of points 
Gi = n n Z" generates the semigroup. The computation of Gi can be done via direct enumeration and 
the knowledge that \Gi\ is the same as the number of cosets of the quotient of Z" by the Abelian group 
generated by the cone generators. 

This way, each simplicial cone ai in the triangulation of C provides us with a set of generators Gi. 
From the union G — LlGi — {wi, . . . ,Wm}, which obviously generates C n Z", we need to find a subset 
H C G whose elements are irreducible and still generate CnZ". The subset H is constructed recursively, 
starting from the empty set, in the fc-th step we check if Wk — h £ C for some h G H. If yes, delete Wk 
from the list and go to the next iteration; otherwise remove all those h in H which satisfy h — Wk 6 C 
and add Wk to H before passing to the next step. Clearly, since we have the inequality representation of 
the cone, it is easy to decide whether a vector belongs to the cone or not. 

With any d-dimensional rational pointed polyhedral cone C = {Ax = 0, a; > 0} and a field k we 
associate a semigroup ring, Rc ^ k[y'^ : a £ Sc\, where there is one monomial y 1^1/2^ ■ ■ ■ y^** in the ring 
for each element a — (oi, . . . , aj) of the semigroup Sc- By the definition of a Hilbert basis we know that 
every element of Sc can be written as a finite linear combination ^ fiihi where the /i^ are nonnegative 
integers. Thus Rc is in fact a finitely generated fc-algebra, with one generator per element of a Hilbert 
bases. Therefore Rc can be written as the quotient k[xi, X2, ■ ■ ■ , xn]/Ic- Once we have the Hilbert basis 
H = {hi, . . . , Hn} for the cone C, Ic is simply the kernel of the polynomial map (j) : k[xi, X2, ■ ■ ■ , xn] — *■ 
k[yi,y2, ■ ■ - yd], where (l){xi) = y''' and for hi = (ai, 02, . . . , ad) we set y'*' = y'l'^yT ■ ■■vT- There are 
standard techniques for computing this kind of kernel (see |^ and references within) . 

It is important to observe that, for our cones of magic arrays, we can give a natural grading to 
Rc- A magic array can be thought of as a monomial on the ring and its degree will be its magic sum. 
For example, all the elements of the Hilbert bases of 3 x 3 magic squares are elements of degree 3. 
Once we have a graded fc-algebra we can talk about its decomposition into the direct sum of its graded 
components Rc = ® -Rc(i)j where each Rc{i) collects all elements of degree i and it is a fc-vector space 
(where i?c(0) = k). The function H{Rc,i) = dimk{Rc{i)) is the Hilbert function of Rc- Similarly one 
can construct the Hilbert- Poincare series of Rc, Hfi^{t) = J2^o ^ : i)t'' ■ 
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Lemma 1.1. Let C be a pointed rational cone, with Hilhert basis H = {hi, . . . , hjsf}. Let the degree of 
a variable Xi in the ring k[xi, . . . ,xn] be the magic sum of the its corresponding Hilbert basis element 
hi. Let Rc be the (graded) semigroup ring obtained from the minimal Hilbert basis of a cone C of magic 
arrays. Then the number of distinct magic arrays of magic constant s equals the value of the Hilbert 
function H{Rc,s). 

Proof: By the definition of a Hilbert basis we have that every magic array in the cone C can be written as a 
hnear integer combination of the elements of the Hilbert basis. The elements of HB(C) = {hi, /i2, . . . , /iat} 
are not affinely independent therefore there are different combinations that produce the same magic array. 
We have some dependencies of the form ^a^/ii = where the sums run over some subsets of 

{1, . . . ,N}. We consider such identities as giving a single magic array. The dependencies are precisely 
the elements of the toric ideal /c, that give Rc = k[xi,X2, ■ • ■ , Xn]/Ic- Every such dependence is a linear 
combination of generators of any Grobner basis of the ideal Jc- Thus, if we encode a magic array X as a 
monomial in variables xi, . . . ,xn whose exponents are the coefficients of the corresponding Hilbert basis 
elements that add to X, we are counting the equivalence classes modulo Ic- These are called standard 
monomials. Finally, it is known that the number of standard monomials of graded degree i equals the 
dimension of Rc{i) as a fc- vector space Chapter 9]. □ 

It is known that the Hilbert-Poincare series of Rc can be expressed as a rational function of the 
form -ffflc {t) = ^l^i_iSi ^ ■ where 5i can be read from the rays of the cone C ; they correspond to the 
denominators of the vertices of the polytope of stochastic arrays whose dilations give the cone C (see 
Theorem 4.6.25 ||2^ and Theorem 2.3 in |2^). To compute the Hilbert-Poincare series we relied on the 
computer algebra package CoCoA that has implementations for different algorithms of Hilbert series 
computations The basic idea comes from the theory of Grobner bases (see §9]). It is known 
that the initial ideal of Ic with respect to any monomial order gives a monomial ideal J and the Hilbert 
functions of k[xi,X2, . . . ,xn]/Ic and k[xi,X2, . . . ,xn]/J are equal. Computing the Hilbert function of 
the monomial ideal J is a combinatorial problem which can be solved by an inclusion-exclusion type 
procedure |^ that eliminates variables at each iteration. 

We illustrate the above algebraic techniques calculating a formula for the number of 3 x 3 magic 
squares, where xq corresponds to the matrix with all entries one, at the bottom of Figure ^, and the other 
4 variables xi,X2,X3, X4 correspond to the magic squares on top of Figure ^, as they appear from left to 
right. The ideal Ic given by the kernel of the map is generated by the two relations xix^—x'^, X2X3—xiXi. 
The first relation means, for example, that the sum of magic square 1 with magic square 4 is the same 
as twice the magic square 5. The CoCoA commands that compute the Hilbert-Poincare series is 

L:=[3,3,3,3,3] ; 

Use S: :=Q[x[l. .5]] ,Weights(L) ; 

I:=Ideal(x[l]*x[4] -x[5] "2 ,x [1] *x [4] -x [2] *x [3] ) ; 
Poincare (S/I) ; 

Non-simplified HilbertPoincare ' Series 

(1 - 2x[l]-6 + x[l]-12) / ( (l-x[l]-3) (l-x[l]-3) (l-x[l]-3) (l-x[l]-3) (l-x[l]-3) ) 

Note that to carry out the computation it is necessary to specify a weight for the variables. In our 
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case the weights are simply the magic sums of the array. It is known that from a rational representation 
like this one can directly recover a quasipolynomial (see pTl §4]). 



|s2 + |s + l if3|s, 

otherwise. 



We have seen already that magic arrays are nonnegative integer solutions of a system Ax = 0,x > 0, 
where A is a matrix with {0, 1, —1} entries. This system defines a pointed rational polyhedral cone C. 
One can set up the cone C as the union of all (real- valued) dilations of the polytope of stochastic magic 
arrays P = {x £ R''' : x > 0, and Bx — 1}. For a positive number n we denote E{P,n) the number of 
lattice points in the dilation nP = {nx\x £ P} = {x E R"^ : x > 0, and Bx = n ■ I}. Note that when n is 
integer, E{P, n) for the P polytope of stochastic magic arrays counts the number of integral magic arrays 
of magic sum n. If we let n take real values, then the union of the different dilations of P as n changes 
is the pointed polyhedral cone C. This is easy to see since any magic square of magic sum A satisfies the 
equations Ax = 0, a: > 0, thus all dilations are contained in the cone C. On the other hand, any solution 
X to the system of inequalities that defines C is a magic square of real valued magic sum A. Dividing all 
entries of the array x by A we obtained a magic array that satisfies the system Bx = 1, x > 0, thus the 
cone C is contained in the union of all dilations of P. It can be verified that the rays of the cone C are 
given by all scalar multiples of vertices of P. For our purposes the main result is a theorem of Ehrhart: 

Lemma 1.2 ( [p!^ ). For a rational k-polytope P embedded in R'', in particular the polytope of stochastic 
magic arrays, the counting function E(P, n) is a quasipolynomial in n whose degree equals k and whose 
period is less than or equal to the least common multiple of the denominators of the vertices of P. 

For example, for 3x3 magic squares the vertices of the polytope of stochastic magic squares are 
obtained by dividing the first 4 magic squares in Figure ^ by 3. In this case the periodicity of the 
function is exactly three. Although in all our computations the period of the quasipolynomial turned out 
to be equal to the least common multiple of the denominators of the vertices of P, this is not true in 
general (see Example 4.6.27 in pTf). 



2 Families of Magic Arrays. 



2.1 Semi-magic Hypercubes: Theorem P7I| 

We consider first the 3x3x3 semi-magic cube. Bona 1^ had already observed that a Hilbert basis must 
contain only elements of magic constant one and two. Here we provide the 12 Hilbert basis elements of 
magic constant 1. There are 54 of magic constant 2, which we are not listing here, but can be downloaded 
from ijww.math.ucdavis . edu/~deloera/RESEARCH/magic .html 



(0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0) 
(1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0) 
(1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0) 
(0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0) 
(1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0) 
(0, 1, 0, 0, 0, 1, 1,0, 0, 1,0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0) 



(0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0) 
(0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1) 
(0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1) 
(0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0) 
(1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1) 
(0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1) 
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From the Hilbert basis and using CoCoA to compute the number of magic cubes we obtain the stated 
rational generating function. Now we claim the number of vertices of the polytope of stochastic semi- 
magic n X n X n cubes is bounded below by (n!)^"/n" . This follows from a bijection between integral 
stochastic semi-magic cubes and n x n latin squares: Each 2-dimensional layer or slice of the integral 
stochastic cubes are permutation matrices (by Birkhoff-Von Neumann theorem), the different slices or 
layers cannot have overlapping entries else that would violate the fact that along a line the sum of the 
entries equals one. Thus make the permutation coming from the first slice be the first row of the latin 
square, the second slice permutation gives the second row of the latin square, etc. From well-known 
bounds for latin squares we obtain the lower bound. 

The polytope of stochastic semi-magic 3x3x3 cubes is actually not equal to the convex hull of 
integral semi-magic cubes. This follows for because the 54 elements of degree two in the Hilbert basis, 
when appropriately normalized, give rational stochastic matrices that are all vertices. In other words, the 
Birkhoff-von Neumann theorem page 108] about stochastic semi-magic matrices is false for 3 x 3 x 3 
stochastic semi-magic cubes. Finally, the polytope of 4-dimensional semi-magic hypercubes has a non- 
integral vertex (each row is a 3-cube worth of values): 



1/3* (0,2, 1,2, 1,0, 1,0, 2, 1,1, 1,0, 1,2, 2, 1,0, 2,0, 1,1, 1,1, 0,2,1 
2,0,1,0,1,2,1,2,0,1,2,0,1,1,1,1,0,2,0,1,2,2,1,0,1,1,1 
1,1,1,1,1,1,1,1,1,1,0,2,2,1,0,0,2,1,1,2,0,0,1,2,2,0,1) 



2.2 Magic Hypercubes: Theorem 

The function that counts magic cubes is a quasipolynomial whose degree is the same as the dimension 
of the cone of magic cubes minus one. For small values (e.g ti = 3,4) we can directly compute this. We 
present an argument for its value for n > 4: 

Lemma 2.1. Let B he the (3n^ + 4) x n'^ matrix with 0, 1 entries determining axial and diagonal sums. 
In this way we see that n x n x n magic cubes of magic sum s are the integer solutions of Bx = 
(s, s, . . . , s)^, a; > 0. For n > 4 the kernel of the matrix B has dimension [n — 1)'^ — 4. 

Proof. It is known that for semi-magic cubes the dimension is {n — Vf which means that the 
rank of the submatrix B' of B without the 4 rows that state diagonal sums is n'^ — (n — 1)'^. It remains 
to be shown that the addition of the 4 sum constraints on the main diagonals to the defining equations 
of the n X n x n semi-magic cube increases the rank of the defining matrix B by exactly 4. 

Let us denote the entries of the cube by . . . ,Xn,n,n and consider the (n— 1) x (rt-- 1) x {n~l) 

sub-cube with entries xi^i^i, . . . , Xn-i,n-i,n-i- For a semi-magic cube we have complete freedom to choose 
these (n—l)^ entries. The remaining entries of the nxnxn magic cube become known via the semi-magic 
cube equations, and all entries together form a semi-magic cube. For example: 

j=l Xi^n,n — 2^i=l 2^j=l ^n,n,n — ^ Z_^i=l 2^j=l 2^k=l ^i,j,k- 
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However, for the magic cube, 4 more conditions have to be satisfied along the main diagonals. Em- 
ploying the above semi-magic cube equations, we can rewrite these 4 equations for the main diagonals 
such that they involve only the variables 2^1,1,1, . . . , Xn-i.n-i.n-i- Thus, as we will see, the complete free- 
dom of choosing values for the variables 2^1,1,1, • ■ • , Xn~i,n-i,n-i is restricted by 4 independent equations. 
Therefore the dimension of the kernel of B is reduced by 4. 

Let us consider the 3 equations in xi,i^i, . . . ,Xn~i,n-i,n-i corresponding to the main diagonals 
Xi^i^ni ■ ■ ■ ,Xn^n,i, Xi^n,i, ■ ■ ■ ,Xn,i,ni and Xn.1.1, ■ ■ ■ ,Xi^n,n- They are linearly independent, since the vari- 
ables Xn-i^n-1,1, 2:„_i.i_„_i, and a;i_„_i^„_i appear in exactly one of these equations. The equation 
corresponding to the diagonal 21,1,1, • ■ • ,Xn,n,n is linearly independent from the other 3, because, when 
rewritten in terms of only variables of the form with 1 < i, j, A: < n, it contains the variable 2:2,2,35 
which for n > 4 does not lie on a main diagonal and is therefore not involved in one of the other 3 
equations. This completes the proof. □ 

We consider now the 3x3x3 magic cubes. There are 19 elements in the Hilbert basis and all of them 
have magic sum value of 3. This already indicates that there is a quasipolynomial counting formula since 
there are no elements of magic sum not divisible by 3. 



(2, 1, 0, 1, 0, 2, 0, 2, 1, 0, 2, 1, 2, 1, 0, 1, 0, 2, 1, 0, 2, 0, 2, 1, 2, 1, 0) 
(1, 2, 0, 1, 0, 2, 1, 1, 1, 1, 0, 2, 2, 1, 0, 0, 2, 1, 1, 1, 1, 0, 2, 1, 2, 0, 1) 
(2, 0, 1, 0, 2, 1, 1, 1, 1, 0, 2, 1, 2, 1, 0, 1, 0, 2, 1, 1, 1, 1, 0, 2, 1, 2, 0) 
(1, 1, 1, 2, 0, 1, 0, 2, 1, 1, 2, 0, 0, 1, 2, 2, 0, 1, 1, 0, 2, 1, 2, 0, 1, 1, 1) 
(1, 2, 0, 2, 0, 1, 0, 1, 2, 2, 0, 1, 0, 1, 2, 1, 2, 0, 0, 1, 2, 1, 2, 0, 2, 0, 1) 
(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1) 
(1, 0, 2, 0, 2, 1, 2, 1, 0, 0, 2, 1, 2, 1, 0, 1, 0, 2, 2, 1, 0, 1, 0, 2, 0, 2, 1) 
(1, 1, 1, 0, 2, 1, 2, 0, 1, 1, 0, 2, 2, 1, 0, 0, 2, 1, 1, 2, 0, 1, 0, 2, 1, 1, 1) 
(1, 0, 2, 1, 2, 0, 1, 1, 1, 1, 2, 0, 0, 1, 2, 2, 0, 1, 1, 1, 1, 2, 0, 1,0, 2, 1) 
(0, 1,2,1, 2, 0, 2, 0, 1, 2, 0, 1, 0, 1, 2, 1, 2, 0, 1, 2, 0, 2, 0, 1, 0, 1, 2) 



(1, 1, 1, 1, 0, 2, 1, 2, 0, 0, 2, 1, 2, 1, 0, 1, 0, 2, 2, 0, 1, 0, 2, 1, 1, 1, 1) 
(2, 1, 0, 1, 1, 1, 0, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 2, 1, 1, 1, 2, 1, 0) 
(2, 1, 0, 0, 2, 1, 1, 0, 2, 1, 0, 2, 2, 1, 0, 0, 2, 1, 0, 2, 1, 1, 0, 2, 2, 1, 0) 
(0, 1, 2, 2, 0, 1, 1, 2, 0, 1, 2, 0, 0, 1, 2, 2, 0, 1, 2, 0, 1, 1, 2, 0, 0, 1, 2) 
(0, 2,1,1, 0, 2, 2, 1, 0, 1, 0, 2, 2, 1, 0, 0, 2, 1, 2, 1, 0, 0, 2, 1, 1, 0, 2) 
(2, 0, 1, 1, 2, 0, 0, 1, 2, 1, 2, 0, 0, 1, 2, 2, 0, 1, 0, 1, 2, 2, 0, 1, 1, 2, 0) 
(0, 2, 1, 2, 0, 1, 1, 1, 1, 2, 0, 1, 0, 1, 2, 1, 2, 0, 1, 1, 1, 1, 2, 0, 1, 0, 2) 
(0, 1, 2, 1, 1, 1, 2, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 0, 1, 1, 1, 0, 1, 2) 
(1, 1, 1, 1, 2, 0, 1, 0, 2, 2, 0, 1, 0, 1, 2, 1, 2, 0, 0, 2, 1, 2, 0, 1, 1, 1, 1) 



From this information, and using CoCoA, we can derive the desired formula for the count that appears 
in Theorem 0.2. Finally we include below an extreme ray for the cone of magic 3x3x3x3 hypercubes. 
Dividing its entries by 15 we get a rational vertex of the polytope of stochastic magic 3x3x3x3 
hypercubes. 



8700877084 475286903 48 10 50 267 
4475286901 10 4852609 10 14285 366 
348 10 50267 10 1 4285366 2 10 3 3 2 10 10 32 



2.3 Magic Squares: Theorem 3.3 



4x4 magic squares: Our calculations using MLP show that there are 20 elements in the Hilbert basis 
for the cone Ca/4x4 of 4 x 4 magic squares. The 8 elements of magic sum one (not 7 as reported in 
and the 12 elements of magic sum 2 are listed below. To save space we present the squares as vectors 

(Xll, . . . , XlA,X21, ■ ■ ■ , X24,, X3I, • • . , X34,, Xu, ■ ■ ■ , Xaa). 



Maya Ahmed et al. 



11 



(0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0) (0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1) (1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0) 
(0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0) (0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0) (1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0) 
(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0) (0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1) 

These 8 permutation matrices are exactly all the magic squares of magic sum 1. The rest of the minimal 
Hilbert basis consists of magic sum 2 magic squares: 



(1,0,1,0,0,0,0,2,0,1,1,0,1,1,0,0) (0,0,2,0,0,1,0,1,1,1,0,0,1,0,0,1) (0,0,1,1,0,1,0,1,1,0,1,0,1,1,0,0) 

(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 2, 1, 0, 1, 0) (1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 2, 0) (0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0) 

(1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 2, 0, 0) (0, 2, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1) (1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1) 

(0, 1, 0, 1, 2, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1) (1, 0, 1,0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1) (0, 0, 1, 1, 0, 1, 1, 0, 2, 0, 0, 0, 0, 1, 0, 1) 



Using CoCoA's Hilbert series computation we obtain the generating function stated in Theorem p.3| . 

5x5 magic squares The 5x5 magic squares are the first challenging case. We were unable so far to 
recover the Hilbert series for this case. By using the fact that the Hilbert basis is a generating set we can 
easily compute several values of the Hilbert function, i.e. the numbers of magic squares for small values 
of the magic sum. Using the generators we consider all possible sums of them with small coefficients, 
making sure that repeated squares are only counted once. The values below allow us to prove that 
there is no polynomial formula that fits those values via interpolation. We use the Ehrhart-Macdonald 
reciprocity laws [ p6[ that give us other 6 values of the function, which together with known roots allow 
for interpolation. There is no solution for the resulting linear system. 



magic sum 


total number magic squares 


1 


20 


2 


449 


3 


6792 


4 


67,063 


5 


484, 419 


6 


2, 750, 715 



The following table lists the number of ele ments in the Hilbert basis. All the elements for all the Hilbert 
bases we have computed can be obtained at rew .math .ucdavis . edu/~deloera/RESEARCH/magic .html 



magic sum 


number of HB elements 


1 


20 


2 


240 


3 


1392 


4 


1584 


5 


1192 


6 


160 


7 


224 


9 


16 




4828 



Finally we prove the rest of Theorem We construct integral extreme ray vectors that, when its 
entries are divided by 1/2, give a fractional vertex of the polytope of stochastic magic squares: Let n > 6 
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and let Pn-2 be an (n — 2) x (rt — 2) permutation matrix that does not contain a non-zero entry on its 
two main diagonals. Let i?„ be the n x n matrix that is constructed as follows: 



• 




= 2 * P„_ 




for i = 2, . . 


• ,n-l, j = 2,.. 


• I ~ li 


• 




— Rn,n,j 


= for i 


= 2,... ,n 


-1, 




• 




— Rn,i,n 


= for i 


= 2,... 


-1, 




• 


Rn,l,l 


— ^n,n,l 






1. 





Since n — 2 > 4, there exists a permutation matrix Pn-2 with no non-zero entries on its main diagonals. 
Thus, Rn is well-defined. 

Lemma 2.2. By construction, Rn is a magic square of size n with magic constant 2, in addition, for 
n > 6, Rn is an extremal ray of the cone of n x n magic squares. 

Proof: Suppose that i?„ is not an extremal ray of the magic square cone. Therefore, there exists a 
non-zero magic square _R„ with magic constant s > whose support is strictly contained in the support 
of Rn- Since every row and column must have at least one non-zero entry, Rn must have a zero in one 
of the corners, that is without loss of generality, Rn.i.i = 0. Since s = X)r=i ^n.i.i = S^=i ^n,n,j, we 

obtain S = Rn,n,l = Rn,n,l + Rn.n.n- ThuS, Rn,n.n — 0. But this COntradictS < S = X]r=l ^n,i,i = 0. 

Therefore, Rn does not exist, implying that i?„ is an extremal ray. □ 



2.4 Pandiagonal Magic Squares: Theorem p74| 

Let us denote by MP„(s) the number of rt x rt pandiagonal magic squares with magic sum s. As in the 
case of magic squares the function MPn{s) is a quasipolynomial in s of degree equal to the dimension 
of the cone plus one. Halleck computed the dimension of the cone to be (n — 2)^ for odd n and 
[n — 2)^ -I- 1 for even n (degree of the quasipolynomial MPn{s) is one less than these). For the 4x4 
pandiagonal magic squares a fast calculation corroborates that there are 8, magic-sum-2, generators. In 
his investigations, Halleck identified a much larger generating set. 



(1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0) (1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0) (0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0) 
(1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1) (0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0) (0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1) 
(1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1) (0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1) 



From the Hilbert basis we can calculate the formula stated in Theorem p.4| using CoCoA. Finally 
we verify that the 5x5 pandiagonal magic squares have indeed a polynomial counting formula. This 
case requires in fact no calculations thanks to earlier work by who proved that for n = 5 the only 
pandiagonal rays are precisely the pandiagonal permutation matrices. It is easy to see that only 10 of 
the 120 permutation matrices of order 5 are pandiagonal: 



(0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0) 
(0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1) 
(0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0) 
(0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0) 
(0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0) 



(1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0) 
(0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0) 
(0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0) 
(1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0) 
(0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1) 
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Once more a simple CoCoA calculation shows that the counting function equals indeed a polynomial, 
g^(s + 4)(s + 3)(s + 2)(s + l)(s^ + 5s + 8){s^ + 5s + 42), as claimed in the Theorem. 
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